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Abstract 



modes that render them useless in numerical integrations of black hole space- 
times, independently of how the equations are differenced. As an example we 
consider the evolution of Schwarzschild and Reissner-Nordstrom spacetimes 
in double null coordinates. 



I. INTRODUCTION 

There are two steps from the Einstein equations to a numerical code for solving them: 
First one selects the fields of the problem (the metric, its derivatives, connection coefficients 
etc.) and a subset of the Einstein equations to evolve these from some boundary conditions. 
This constitutes the "differential problem". The resulting set of differential equations is 
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then transformed into a set of difference equations on a numerical grid, constituting the 
"difference problem". (Here the term "boundaries" includes spacelike and null as well as 
timelike boundaries.) 

In the first step, one distinguishes between "free" and "constrained" evolution schemes. 
In a free evolution scheme, the constraints are imposed only at the boundaries, and all fields 
are propagated by means of the evolution equations. In a constrained scheme, one uses 
fewer evolution equations, and instead reconstructs some fields at each time step from the 
constraints. (Here the term "time" refers to the coordinate that labels slices in the numerical 
evolution, which can be null. Constraints are all equations restricted to a hypersurface.) 

In this paper we show that free evolution may already be ill-posed as a differential 
problem in some cases because it admits exponentially growing perturbations which are 
solutions of the evolution equations, but which violate the constraints and are therefore 
unphysical. For analytic purposes this does not matter, as imposing the constraints on 
the boundaries assures that they are obeyed everywhere, but any finite-differencing scheme 
introduces small perturbations which do not generically obey the constraints, and serve 
as initial data for the growing unphysical modes. The amplitude of these modes depends 
therefore on the choice of difference scheme (and decreases with the grid spacing), but their 
evolution (for example, the rate of exponential blow-up) is determined by the differential 
problem alone. 

As an example differential problem we consider free evolution in spherical symmetry in 
double null coordinates. It is sufficiently simple that we can calculate the blowup of the 
unstable mode analytically, before confirming it in a testbed calculation. This particular 
instability is a serious problem inside and just outside black holes, harmless far outside black 
holes, and is not present in perturbations around flat spacetime. 

II. DOUBLE NULL COORDINATES IN SPHERICAL SYMMETRY 

The metric of any spherically symmetric spacetime can be written as 
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ds 2 = -Afdudv + r 2 dn 2 , (1) 

where the metric coefficients / and r depend on u and v only. Our matter will be a minimally 
coupled massless real scalar field 0, plus a spherically symmetric Maxwell field. The Maxwell 
field has no sources, and is therefore constrained to be a pure Coulomb field of constant 
charge q anchored at (the singularity) r = 0. 

This model has been used in both analytic ]IJ and numerical |2|,|3|] work as a toy model for 
studying the interior of realistic black holes. Realistic black holes should be spinning, at least 
slightly, and would therefore be expected to have a Cauchy horizon, like the Kerr solution. 
On the other hand they would be inevitably perturbed by gravitational wave tails, and these 
perturbations would blow up on the Cauchy horizon. What really happens has therefore 
been the subject of prolonged investigation. Replacing the rotation by an electric charge 
and gravitational waves by the scalar field allows the simplification of spherical symmetry. 
Although we have chosen a particular matter model here, we shall see that the instability 
arises in the gravitational part of the equations and is therefore independent of the matter 
model. 

The uu and vv components of the Einstein equations are 

r,„„-^r, tt + (^„) 2 = ^i = 0, (2) 
r, vv ~ y- r, v + (<fr :V ) 2 = E 2 = 0. (3) 



The 99 component is 



. r ,u r ,v f . Q f u i\ I i \ 

r, uv H 1 5- = Hi = 0, (4) 



and the uv component is 



+ 2^ - 2^ + 2<P^ V = H 2 = 0. (5) 

The other components of the Einstein equations are redundant. The massless wave equation 
for 0, restricted to spherical symmetry, is 



<P,uv + —<t>,u + —<t>,v = H 3 



= 0. 



(6) 



We then have three hyperbolic equations Hi, H 2 and H 3 and two elliptic equations Ei and 
E 2 . (The elliptic equations are really ordinary differential equations because of the spherical 
symmetry.) The elliptic equations are propagated by the hyperbolic equations in the sense 



A possible evolution scheme, perhaps the most natural one, and certainly one easy to 
implement numerically, is to consider the three hyperbolic equations as evolution equations 
for /, r and <fi, and the elliptic equations as constraints which are imposed only on the 
boundary. The natural initial value problem for these equations, is a double null initial 
value problem, r, / and <fi are given as functions of u on the null cone v = v , subject to the 
constraint E\ = 0, and as functions of v on the intersecting null cone u = uq, subject to the 
constraint E 2 = 0. At the intersection point (u — uq,v — vq) it is sufficient that the data 
r, / and <fi be continuous. In the free evolution scheme, the constraints are imposed only 
on the null boundary data, but are not used for the numerical solution inside the numerical 
domain. 

We now examine the stability of this free evolution scheme in the context of black hole 
physics. As a testbed case we use the Reissner- Nordstrom solution. It is the unique solution 
for = 0, and it is known in double null coordinates in (essentially) closed form. In one 
particular gauge choice (Eddington-Finkelstein coordinates) this solution is 



of 



Ei, v = H hu + 2r0 iU #3 - r -^Ei - r, u H 2 + 



(7) 



while a similar equation holds for E 2 . 



III. THE FREE EVOLUTION SCHEME 




(8) 



where r{u,v) is given implicitly by 
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r(u,v)=r(r*), r* = ^7 = f(r). (9) 

The implicit equation dr/dr* = f(r) can be solved numerically to arbitrary precision so that 
the solution exists in closed form for all numerical purposes. 

In our test case we have set = because we expect a nonvanishing to give rise to a 
physical instability, namely mass inflation, inside the black hole, while we want to check if 
such instabilities are already contained in the free evolution scheme in a situation when we 
know that no physical instabilities can be present. 

In order to look for instabilities in analytic approximation, it is useful to make a change 
of variables in order to eliminate first derivatives in the hyperbolic equations. With the new 
variables 

y = m/ + mr, x = r 2 (10) 

the evolution equations Hi — and H 2 = 0, now restricted to = 0, become 

x, uv + A(x,y)=0, y,uv + B(x,y) = 0. (11) 

(The evolution equation H 3 = for is dropped.) The constraints E\ — and E 2 = 
become 

uu y,u%,u 0, X t vv y tV x jV 0. (12) 

Now we linearize the evolution equations, denoting the perturbations of x and y by £ and r\: 

C,uv + A tX £ + A >y r) = 0, r) >uv + + B !y r) = 0. (13) 

The coefficients of these equations on the Reissner-Nordstrom background are 

a = 4„ = -2/(i-£), 5 = ^ = ^ = ^(1-3^, ^ = -^(1-5^). 

(14) 

To obtain an analytic approximation to the perturbations we now consider r and / of 
the background as slowly varying functions of u and v. We make a mode ansatz 



The ansatz turns the derivative £ )nt , into the algebraic expression (k 2 — uj 2 )^k, and the 
equations (|T3|) into the local dispersion relation oj 2 {k) = k 2 + A±. The X± are the eigenvalues 
of the two-by-two matrix [(A jX , A^ y ), (B >x , B jV )]. Their value on the Reissner-Nordstrom 
background is 



/ 



^± — 



1 - 3p)±V3(l -/£>)(! -5p) 



where p = q jr , (16) 



and the corresponding eigenvectors are 



7]k± = ± 2^i i—p-^- (17) 
With u and v both increasing to the future, u + v labels time, and u — v labels space. 
u + v increases away from the null boundary data. (Inside the black hole, evolution towards 
increasing u + v is also evolution towards the singularity.) Therefore we have an instability 
if cu 2 (k) is negative for any k, that is, for A < 0. We see that for q — 0, A± = (1 ± \^3)r~ 2 f. 
Although / changes sign at the horizons, one of the A is always negative for all r. For q ^ 
the situation is not changed drastically. There is now a region where both A are positive, 
namely 0.74 \q\ < r < \q\. But this is only a narrow band inside the black hole, and globally 
the free evolution scheme is still unstable. 

The transformations that leave the form ([!]) of the metric invariant are u — > U(u) and 
v —>■ V(v), where U and V are arbitrary functions. The speed of the exponential blowup 
itself is gauge- invariant, in the sense that A transforms correctly under the coordinate trans- 
formations u — > U(u) and v — > V(v). We can see this noting, from ([I]), that / transforms in 
the same way as g jUV for any scalar g, and that A is the product of / times a scalar (r is a 
scalar), and hence also transforms like /. Our results would therefore also hold in Kruskal 
coordinates, or any other double null coordinates. 

We should stress again that the unstable perturbations we have constructed are not tied 
to a particular numerical scheme. They are solutions to the hyperbolic part of the Ein- 
stein equations. They must arise in any numerical algorithm implementing a free evolution 
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scheme in null coordinates, simply because they are excited by the discretisation error of 
any numerical scheme. 

Of course these unstable modes are not solutions of the constraints, or elliptic part of 
the Einstein equations. We can explicitly verify this by constructing the perturbations that 
do obey the constraints. The linearized constraints are 

Because the Reissner-Nordstrom solution is unique, any perturbations of (HH) obeying all 
the Einstein equations must be infinitesimal coordinate transformations. Perturbatively we 
write u — > u + fi(u) and v — > v + u(v). The corresponding perturbative changes in the metric 
variables are 

£ = -x )U ft - x ;V v, r) = -y iU fi-y iV v-dfi/du-dv/dv, (19) 

Clearly these obey the linear constraints and are bounded. The unstable modes we have 
given do not obey the linear constraints and blow up. 

We first encountered this instability in a code modeled on that of reference |2|, when 
we found we were unable to recover the Reissner-Nordstrom solution in a testbed. Having 
constructed the unstable modes analytically, we could quantitatively verify their presence in 
the code. (Once more it should be said that the instability is connected to the free evolution 
scheme, not to any particular numerical implementation.) As an example we numerically 
construct a null diamond (that is, a square in u and v) of the Schwarzschild solution centered 
on r = 3M. (We have chosen r ~ 3M because there the instability is greatest outside 
the horizon.) Fig. 1 shows the numerical setup. At r = 3M in the coordinates (§,0), 
A + ~ — 0.025 M -2 . This corresponds to a blowup of the unstable modes as exp(0.32 M _1 t), 
where t — (u + v)/2 is the usual Schwarzschild time coordinate. Numerically we find that 
the error grows as exp(0.24 M~ x t) at r = 3M. The discrepancy of the exponent may be due 
to the approximation of constant r and / we have made in the analytic calculation. 

There is a second prediction of the analytic model which can be verified at least qualita- 
tively. From the form of the eigenvectors fll7D , together with the definitions (|T0|) it follows 
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that the relative error in / is generically of the same order as the relative error in r, except 
for r ~ \q\, where it will much greater (because 1 — p <C 1). This feature is confirmed by 
numerical evolutions inside and outside the horizon, for various values of q. 

IV. CONCLUSIONS 

We have given a simple clean example of a numerical instability arising from constraint 
violation. The amplitude of the unstable mode depends on the discretisation scheme, but its 
growth rate does not. We have estimated the growth rate analytically, and have confirmed 
it for a particular discretisation. 

The free evolution scheme we have described was used for the numerical study of per- 
turbed Reissner-Nordstrom black holes in reference 0. As our investigation shows, this 
scheme is unable even to evolve the exact Reissner-Nordstrom solution (setting the scalar 
field to zero) because of an exponential instability eventually crashing the code. It can 
therefore not be used for the study of the physical instability triggered by small physical 
perturbations. (This is the physical instability [Q] predicted to lead to mass inflation and the 
destruction of the Cauchy horizon.) This casts a shadow on the physical results obtained 
using this scheme. 

Other published codes using a double null numerical grid use fully constrained 

evolution. The interior of a charged black hole has been treated with a fully constrained 
evolution scheme based on coordinates u (retarded time) and r (curvature radius), but 
evolving them on a double null grid ||. This kind of algorithm has been pioneered in [|J, 
and has been used extensively since 0-0. A code actually based on double null coordinates, 
as well as a null grid, was used in ||, implementing a fully constrained scheme. Here 
fully constrained means that one uses the maximum number of equations containing only 
^-derivatives. 

Another instance of an ill-posedness and instability arising from constraint violation has 
been found for the nonlinear evolution of weak gravitational waves on a flat background || . 
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The growth rate of the instability depended on the discretisation scheme, however, which is 
puzzling in the light of our results. 

Generally our results suggest that free evolution schemes are generally harder (or impos- 
sible, as in this case) to make stable than fully constrained schemes. This is not in conflict 
with the statement [TtJ that if one if one evolves in free evolution to a certain order in 
the grid spacing, the constraints are automatically obeyed to that order. The numerical, 
constraint-violating, error may go down as some power of decreasing step-size, but it can 
also be growing exponentially with time. If this exponential growth is fast enough, one will, 
in typical situations, be unable in practice to compensate for the error at late times by using 
a finer numerical grid. 

Enforcing all the constraints at each time step can avoid a blowup only when the solu- 
tion to be calculated is itself insensitive to small perturbations in the initial data, because 
then all rapidly growing perturbations must be constraint violations. The solution itself 
however may be highly sensitive to perturbations in the initial data, and/or may blow up 
at a spacetime singularity. Examples of such physical instabilities are mass inflation in the 
perturbed black holes we discussed here, the threshold of black hole formation, or chaos 
in the mixmaster universe [ |TT] , P~2] | . Such examples may not be strictly well-posed but are 
nevertheless interesting. Then one must enforce the constraints to distinguish physical from 
unphysical (constraint-violating) perturbations. 
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1 ^ 

r=3m 

FIG. 1. The numerical grid setup for our example. The partially filled-in numerical grid indi- 
cates in what order grid points are completed. The fat lines u = uq and v = vq are null boundaries. 
From a numerical point of view one might call the former the initial time and the latter a bound- 
ary, but clearly this is a matter of words. The dotted line marks the boundary of the domain of 
dependence of the null boundary data, t and r* are the Schwarzschild time and tortoise radial 
coordinate. 
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